The nonequilibrium discrete nonlinear Schrodinger equation 
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We study nonequilibrium steady states of the one-dimensional discrete nonlinear Schrodinger 
equation. This system can be regarded as a minimal model for stationary transport of bosonic 
particles like photons in layered media or cold atoms in deep optical traps. Due to the presence of two 
conserved quantities, energy and norm (or number of particles), the model displays coupled transport 
in the sense of linear irreversible thermodynamics. Monte Carlo thermostats are implemented to 
impose a given temperature and chemical potential at the chain ends. As a result, we find that the 
Onsager coefficients are finite in the thermodynamic limit, i.e. transport is normal. Depending on 
the position in the parameter space, the "Seebeck coefficient" may be either positive or negative. For 
large differences between the thermostat parameters, density and temperature profiles may display 
an unusual nonmonotonic shape. This is due to the strong dependence of the Onsager coefficients 
on the state variables. 

PACS numbers: 05.60.-k 05.70.Ln 44.10.-K 



I. INTRODUCTION 



The Discrete Nonlinear Schrodinger (DNLS) equation 
PI has important applications in many domains of 
physics. As it is well known, such equation arises in sev- 
eral different problems. A classical example is electronic 
transport in biomolecules [|| . In the context of optics or 
acoustics it describes the propagation of nonlinear waves 
in a layered photonic or phononic system. Indeed, in 
a suitable limit, the dynamics of high-frequency Bloch 
waves is described by a DNLS equation for their enve- 
lope (see Refs. 0, Q for details). On the other hand, in 
the realm of the physics of cold atomic gases, the equa- 
tion is an approximate semiclassical description of bosons 
trapped in periodic optical lattices (see e.g. Ref. Q and 
references therein for a recent survey) . Many other phys- 
ical problems have been recently addressed having the 
DNLS equation as a basic reference model, like the ef- 
fect of nonlinearity on Anderson localization 7, Q and 
the violation of reciprocity in wave scattering 9( just to 
mention a few recent examples. 

While a vast literature has been devoted to localization 
problems, much less is known about finite-temperature 
properties. The first analysis of the equilibrium statis- 
tical mechanics of DNLS systems has been performed 
in Ref. pi| , while the relaxation of localized modes (dis- 
crete breathers) in the presence of phonon baths has been 
discussed in [TJ [l2j]. Several results can be translated 
to other types of nonlinear lattices, where a DNLS-likc 
equation represents an approximation of the lattice dy- 



namics (l 31 ] . 

An even less explored field is that of nonequilibrium 
properties of the DNLS equation In particular, the 
case of an open system that exchanges energy with exter- 
nal reservoirs has not been treated so far. The presence 
of two conserved quantitities naturally requires to argue 
about coupled transport, in the sense of ordinary lin- 
ear irreversible thermodynamics. Despite the very many 
studies of heat conduction in oscillator chains [H, fla j. 
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works on coupled transport are scarce [17H19I ] . Interest 
in this field has been revived by recent works on ther- 
moelectric phenomena [13, [H[ in the hope of identifying 
dynamical mechanisms that could enhance the efficiency 
of thermoelectric energy conversion [22|, [23[ . 

In order to investigate transport properties, we need 
to introduce the interaction of the system with external 
reservoirs that are capable to exchange energy and/or 
norm. For models like DNLS this is much less straight- 
forward than for standard oscillator chains, where e.g. 
Langevin thermostats are a typical choice plj . Here 
we propose and test a very simple Monte Carlo scheme 
which is easy to implement and suitable for the model at 
hand. Another important difference between the DNLS 
and standard oscillator chains (like the Fermi-Pasta- 
Ulam or Klein-Gordon models) is that its Hamiltonian 
is not the sum of kinetic and potential energies. Thus, 
it is necessary to introduce suitable operative definitions 
of kinetic temperature T and chemical potential fj,, to 
measure such quantities in actual simulations. In the fol- 
lowing, we make use of a recent definition of the micro- 
canonical temperature [24j and extend it for the estimate 
of the chemical potential. 

By imposing small T and fi jumps across the chain, we 
can determine the Onsager coefficients, which turn out 
to be finite in the thermodynamic limit, i.e. both en- 
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ergy and mass conductions are normal processes. From 
the Onsager coefficients we can thereby determine the 
"Seebeck coefficient" S [25| which we find to be either 
positive or negative, depending on the thermodynamic 
parameters (i.e., mass and energy density). For larger 
temperature or chemical-potential differences, although 
one can still invoke the linear response theory, some sur- 
prising phenomena emerge. One example is the "anoma- 
lous heating" that can be observed when the chain is 
attached to two thermostats operating at the same tem- 
perature: along the chain, T reaches values that are even 
three times larger than that imposed on the boundaries. 
This phenomenon can be observed only in the case of 
coupled transport, since it is due to the variable weight 
of the non-diagonal terms of the Onsager matrix. It is 
apparent in the DNLS, because of the strong variability 
of the Onsager coefficients. 

The paper is organized as follows. In Sec. [II] we intro- 
duce the model and describe the heat baths. In Sec. IIII1 
we define the relevant thermodynamic observables and 
the formalism (e.g., the Onsager coefficients) necessary 
to characterize nonequilibrium steady states. Sec. IIVI is 
devoted to a discussion of the stead states, both in the 
case of small and large T, fi differences. In Sec. [V] we 
provide a pictorial representation of the general trans- 
port properties, by reconstructing the zero-flux curves. 
Finally, the last section is devoted to the conclusions and 
to a brief summary of the open problems. 



the grand-canonical ensemble with the help of trans- 
fer integral techniques. It is schematically described in 
Fig. [TJ the lower dashed line corresponds to the ground 
state (T = 0) upon varying the particle density; the 
upper dashed line corresponds to infinite temperature. 
The nonequilibrium studies described in this paper cor- 
respond to the region in between such two curves. 
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FIG. 1: (Color online) Parametric plots of the local norm and 
energies [a(y), h(y)] for T L = T R = 1, fi L = 0, (ir = 2 and 
increasing chain lengths N = 200, 800, 3200 (solid lines, bot- 
tom to top). The three (blue) dashed lines are the isothermal 
T = 0, T = 1 and T = oo respectively. Lines at constant 
chemical potential (open symbols) y, = 0, fi = 1 and /i = 2 
(left to right respectively) are obtained by equilibrium simu- 
lations. 



II. SETUP 

In one dimension, the DNLS Hamiltonian writes 

^ N N-l 

i=l i=l 

where the sum runs over the N sites of the chain. The 
sign of quartic term is positive, as we refer to a repulsive- 
atom BEC, while the sign of the hopping term is irrel- 
evant, due to the symmetry associated with the canon- 
ical (gauge) transformation z n — > z n e™ n (where z n = 
(p n + iq n ) /y/2 denotes the amplitude of the wave func- 
tion). The equations of motion are 

^n+1 Zn — 1 2"|z n | Z n (2) 

with n = 1, • • • ,N, and fixed boundary conditions (zq — 
zn+i = 0). The model has two conserved quantities, the 
energy and the total norm (or total number of particles) 

JV 

A = 5>? + a . (3) 

8=1 

As a consequence, the equilibrium phase-diagram is two- 
dimensional, as it involves the energy density h = H/N 
and the particle density a = A/N. The first reconstruc- 
tion of the diagram was carried out in Ref. [lOf within 



We aim at characterizing the steady states of the chain 
when put in contact (on the left and right boundaries) 
with two thermostats at temperature Tl and Tr and 
chemical potentials hl an d /xr, respectively. The im- 
plemcntaton of the interactions with a heat bath is of- 
ten based on heuristics. In particle models, the sim- 
pler schemes consist in either adding a Langevin noise, 
or in assuming random collisions with an equilibrium 
gas [HI, [lH . For the DNLS this is less straightforward: 
adding white noise and a linear dissipation drives the 
system to infinite temperature, i.e. to a state in which 
relative phases are uncorrelated. 

In the absence of a first-principle definition of heat 
bath, we consider two phenomenological Monte- Carlo 
heat baths. The general scheme of this kind of heat bath 
involves a stochastic dynamics which perturbs the canon- 
ical variables p\ — > p\ + Spi and q\ — > qi + 8q\ [26[ at 
random times, chosen according to a uniform distribution 
in the interval [t m i n , t max \. The perturbations Sp and Sq 
are independent random variables uniformly distributed 
in the interval [-R, R]. Moves are accepted according to 
the standard Metropolis algorithm, evaluating the cost 
function exp {-T^Aif - fi L AA)} with T L and [i L be- 
ing the temperature and the chemical potential of the 
heat bath. This kind of thermostat exchanges both en- 
ergy and particles. In some cases, however, we need to 
study the chain behavior in the absence of one of the two 
fluxes (energy and norm). A simple way to study these 
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setups is to modify the perturbation rule of the thermo- 
stat, requiring the exact conservation of the correspond- 
ing local density (energy density or norm density). We 
have thus the following two schemes: 

Norm conserving thermostat- The perturbation acts 
only on the phase (f>i of the complex variable Z\. More 
precisely we impose 4>\ — > <f>i +S(f>i mod{2'K), where 5<p is 
a random variable, uniformly distributed in the interval 
[0,27r]. This dynamics conserves exactly the local ampli- 
tude \zi\ 2 and therefore the total norm A. 

Energy conserving thermostat- In this case, it is neces- 
sary to go through two steps to conserve the energy 



where ( ) stands for the microcanonical average, 



ei = |zi| 4 + 2|^i||2 2 | cos ( 



(4) 



First, the amplitude \z\\ is randomly perturbed. As a 
result, both the local amplitude and the local energy 
change. Then, by inverting, Eq. (0}, a value of 4>i that 
restores the initial energy is seeked. If no such solution 
exists, we go back to the first step and choose a new 
perturbation for | z\ \ . 

There is a basic difference between the two types of 
thermostats. In the general scheme, a steady state is 
characterized by four parameters Tl, Tr, hr. On 
the other hand, for the norm-conserving scheme we only 
assign Tl, Tr and the norm density atot of the whole 
chain. As a consequence, the value of /x on the bound- 
ary is not fixed and must be computed from the simu- 
lation. If the steady state is unique, the former ther- 
mostating scheme must yield the same results, once the 
chemical potentials are suitably fixed. A numerical test 
of this equivalence has been performed, by reconstruct- 
ing some zero-flux profiles with both thermostats. The 
curves overlap reasonably well, although some small sys- 
tematic deviations are present. This is because the norm 
flux is never exactly zero in the non-conservative case 
(typically of order ~ 10~ 4 in a chain of 1000 sites). In 
addition, there are slightly different thermal resistance 
effects in the two schemes. Besides those discrepancies, 
we conclude that the proposed schemes work equally well 
for the generation of nonequilibrium steady states. 



III. PHYSICAL OBSERVABLES 

In order to characterize the thermodynamic properties 
of the DNLS, we extend the approach of Ref. [24| to de- 
rive an operative definition not only of the microcanoni- 
cal temperature but also of the chemical potential. The 
starting point are the usual definitions T -1 = dS/dH, 
and fi/T = —dS/dA, where iS is the thermodynamic en- 
tropy. The partial derivatives must be computed taking 
into account the existence of two conserved quantities 
(hereafter called C\ and C 2 ). Thus, 



as 



U\\w y 



(5) 



w 2 



VCi (VCi • VC 2 )VC 2 
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2N 

E 
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IVdllHVCa 



(6) 



dCi dC 2 dCi dC 2 
dxj dxk dxk dxj 



and X2j = Qj, x 2 j+i = Pj. By setting C\ — H and 
C*2 = A, the above formula reduces to the expression for 
T derived in [24|. Moreover, by assuming C\ — A and 
C 2 = H, Eq. ([5]) defines the chemical potential fi. Notice 
that both expressions are nonlocal. Nevertheless, we have 
verified that it is sufficient to compute the expression 
§5§ over as few as 10 sites to obtain, after some time 
averaging, reliable "local" estimates of both T and \i [13] • 
The expressions for the local energy- and particle- 
fluxes are derived in the usual way from the continuity 
equations for norm and energy densities, respectively 



3a(n) 
jh(n) 



2 (Pn+iqn - 
- (PnPn-1 



PnQn 



-i) 

-i) 



(7) 

(8) 



The approach to the steady state is controlled by veri- 
fying that the (time) average fluxes are constant along 
the chain (j a (n) = j a and jh{n) = jn)- Moreover it is 
also checked that j a and jh are respectively equal to the 
average energy and norm exchanged per unit time with 
the reservoirs. 

As usual in nonequilibrium molecular dynamics simu- 
lations, some tuning of the bath parameters is required 
to minimize boundary resistance and decrease the statis- 
tical errors, as well as the finite-size effects [l5j]. For our 
Monte-Carlo thermostats, we observed that it is neces- 
sary to tune the perturbation amplitude R. Typically, 
there is an optimal value of R for which one of the two 
currents is maximal (keeping the other parameters fixed) , 
but this value may depend on T and (i. Since it would 
be unpractical to tune the thermostat parameters in each 
simulation, we decided to fix them in most of the cases. 
In particular we have chosen R = 0.5, t m i n = 10~ 2 and 
tmax = 10 _1 . Some adjustments have been made only 
when the fluxes were very small. 

In the thermodynamic limit (i.e. for sufficiently long 
chains), the local forces acting on the system are very 
weak and one can thereby invoke the linear response the- 
ory. This means that forces and fluxes are related by the 
relations 1231 



Ja ^aa , "T J^ah , 

ay ay 



(9) 



jh — —Lha 



dim 

dy 



L 



hh 



dp 
dy 



where we have introduced the continuous variable y — 
i/N, while /3 denotes the inverse temperature 1/T; L 
is the symmetric, positive definite, 2x2 Onsager ma- 
trix. Notice that the first term in the r.h.s. of the above 



equations is negative, since the thermodynamic forces are 
— /3/i and li and that detL = LaaL^ — L\ a > 0. 

The particle (a) and thermal (k) conductivity can be 
expressed expressed in terms of L, 



> detL 



(10) 



Analogously, the Seebeck coefficient S, which corre- 
sponds to (minus) the ratio between the chemical- 
potential gradient and the temperature gradient (in the 
absence of a mass flux), writes 



,s'= ;('^ 

J-'n.n. 



/' 



(11) 



We conclude this Section, by mentioning another im- 
portant parameter, the figure of merit 



ZT 



aS 2 T (L ha - fiL a 



detL 



which determines the efficiency rj for the conversion of a 
heat current into a particle current as [23[ 



V = Vc 



^JZT +1-1 



y/ZT +1+1 



For large ZT, rj approaches the Carnot limit rjc- Un- 
derstanding the microscopic mechanisms leading to an 
increase of ZT is currently an active topic of research 



IV. STEADY STATES 
A. Local analysis 

In a first series of simulations we have studied the 
nonequilibrium states in the case of small differences be- 
tween the two thermostats, verifying that transport is 
normal, i.e. the Onsager coefficients are finite in the 
thermodynamic limit. This is less obvious than one could 
have imagined |28[ . In any case, for fixed AT = Tr — T L 
and A/i = /iR — /it, the two fluxes j a and jh are inversely 
proportional to the system size N. At high enough tem- 
peratures, the asymptotic scaling sets in already in chains 
a few hundred sites long (see Fig. [5^,). Moreover, if AT 
A/it and are small enough, the profiles of T and ji along 
the chain are linear as expected. 

However, upon decreasing the temperature, the mini- 
mal chain length needed to observe a normal transport, 
becomes very large. As shown in Fig. [2j3, for the same 
range of lattice sizes as in panel a, the currents are al- 
most independent on TV, as one would expect in the case 
of ballistic transport. This is because at small temper- 
atures, one can always linearize the equations of motion 
around the ground state (which depends on the norm 
density), obtaining a harmonic description and thereby 
an integrable dynamics. 




FIG. 2: Average energy current (squares) and norm current 
(dots) versus chain size N: (a) High-temperature regime Tl = 
2, Tr — 4, /i = and (b) Low-temperature regime Tl = 0.3, 
T R = 0.7, ji = 1.5 The thin (blue) line is the 1/N behavior 
expected for normal transport. Each value is obtained by 
computing the fluxes on a run of 5 x 10 6 time units. 



A plot of the four Onsager coefficients in the (T, fx) 
plane is reported in Fig. [3J Within statistical errors, the 
off-diagonal terms are always positive in the considered 
range. All coefficients are larger for small T and large li. 
This is connected to the scaling behaviour of the linear 
coefficients in the vicinity of the ground state [28{ . 

The resulting coefficient S is plotted in Fig. where 
one can see that there are two regions where the See- 
beck coefficient is positive, resp. negative, separated 
by a curve which, according to Eq. (fTTI) . is defined by 
Lha/Laa = M (see below). This means that the relative 
sign of the temperature and chemical-potential gradients 
is opposite in the two regions (in the presence of a zero 
norm-flux) . This is indeed seen in Fig. \%jp where the re- 
sult of two different simulations are plotted in the two 
regions. 

Finally, since the figure of merit ZT roughly follows 
S, there is only a modest change in the considered pa- 
rameter ranges. Moreover, for fixed T, ZT decreases 
upon increasing ji. This is qualitatively in agreement 
with the general expectation that an increasing strength 
of interaction (increasing ji means increasing the average 
norm and thus the nonlinearity) is detrimental for the 
efficiency. 



B. Global analysis 

If the temperature- or the chemical-potential differ- 
ence is no longer small, the temperature and chemical- 
potential profiles are expected to have a nonlinear shape. 
This is because, as we have seen in the previous subsec- 
tion, the Onsager matrix varies with a and h (or, equiv- 
alently, with T and ll). 
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FIG. 3: Elements of the Onsager matrix L in the T, /i plane, 
for a chain of length N = 500; AT = 0.1, Afi = 0.05. Each 
value is obtained by computing the fluxes on a run of 5 x 10 6 
time units. 




A particularly striking example is reported in Fig. [5] 
Both T(y) and n(y) profiles do approach the imposed 
values at the chain edges (up to tiny jumps due to the 
boundary impedance). However, T(y) exhibits a remark- 
able non-monotonous profile: although the chain is at- 
tached to two heat baths with the same temperature, it 
is substantially hotter in the middle (up to a factor 3!). 

Another way to represent the data is by plotting the 
local norm and energy densities in the phase plane (a, h). 
By comparing the results for different chain lengths, we 
see in Fig. Q] that the paths are progressively "pushed" 
away from the T = 1 isothermal and for TV = 3200 the 
asymptotic regime is attained. 

In order to understand the onset of such anomalous 
shape, it is convenient to rewrite Eq. (|9]) by referring to 
T and fi. By introducing vector notations, we can write, 



J = A( M ,T)^ 
ay 



(12) 



where J = {j a ,jh), v = (/i, T), while the matrix A (which 
is no longer symmetric) can be expressed in terms of the 
Onsager matrix and of the fields T and \x (for instance, 
An = —L aa /T). By now inverting the above equation 
one obtains 



dv 

dy 



= A- X ( M ,T)J 



(13) 



where A -1 denotes the inverse of A. This system de- 
scribes a set of two linear differential equations which are 
non-autonomous (since the matrix coefficients in general 
vary with /j, and T). 



FIG. 4: (Color online) (a) Seebeck coefficient S obtained from 
the data in Fig. [3] (b) Temperature and chemical proten- 
tial profiles for Tl = 1, Tr = 1.5; simulation with norm- 
conserving thermostats at two values of the norm density atot 
corresponding to values of S with opposite signs. 




FIG. 5: (Color online) (a) Temperature and chemical poten- 
tial profiles as a function of y = i/N for a chain of iV = 3200 
sites and Tl = Tr = 1, = 0, (ir = 2. Each point is an 
average of the appropriate microcanonical expression derived 
from Eq. ([5]) over a subchain of about 10 sites around i. (b) 
Norm and energy densities corresponding to the profiles in 
(a). 
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If one assumes to know the "material" properties (i.e. 
the matrix A~ x ) and wishes to determine fluxes and pro- 
files, can proceed by integrating the differential equa- 
tions, starting from the initial condition T(0) = Tl, 
/j,(0) = lil- The, a priori unknown, parameters j a and jh 
can be determined by imposing that the final condition 
is T(l) = Tr and fi(l) = lir. Alternatively, if the fluxes 
are known, one can integrate the equations up to any 
point 2/0 ■ and thereby generate the profiles that would be 
obtained by attaching the right end of the chain to ther- 
mal baths with temperature Tr — T(tjq) and chemical 
potential ll r = ii(y ). 

In order to check the validity of the method, we have 
also adopted an alternative point of view, by combin- 
ing the knowledge of the fluxes with simulations of short 
chains and small gradients to determine the elements of 
the matrix A in suitably selected points in the (T, /i) 
plane. In order to estimate the four entries of A, it is 
necessary to perform two independent simulations for, 



T LtR — T ± AT 



Tl,r = T 
li l ,r = n ± A/x 



With such information, we have been able to estimate 
dvi/dy along the chain (from Eq. (|13p ) and to compare 
the results with the direct simulations. The results plot- 
ted in Fig. [6] demostrate that the two approaches are in 
excellent agreement. 




FIG. 6: (Color online) Spatial derivatives of T and \x (panels 
(a) and (b), respectively) computed from the profiles of Fig. [5] 
(black dashed lines) and their reconstruction (red solid lines) 
by Eq. (|13[) using the linear response coefficients (matrix A). 
The latter have been calculated on a chain of N — 250. The 
quality of the recostruction improves by increasing the lattice 
size has shown by tha blue filled dot which is obtained for 
N = fOOO. 



vanishing fluxes j a and jh- They can be directly de- 
termined by means of the conservative thermostats pre- 
sented in Sec. II. Some lines are plotted in Fig. [7] both in 
the plane (a, h) and (T, li). It is worth recall that in the 
absence of a mutual coupling between the two transport 
processes (zero off-diagonal elements of the Onsager ma- 
trix) such curves would be vertical and horizontal lines 
in the latter representation. It is instead remarkable to 
notice that the solid lines, which correspond to jh = 
are almost veritical for large li: this means that in spite 
of a large temperature difference, the energy flux is very 
small. This is an indirect but strong evidence that the 
nondiagonal terms are far from negligible. 

The condition of a vanishing particle flux j a = defines 
the Seebeck coefficient which is S = —dfi/dT. Accord- 
ingly, the points where the dashed curves are vertical in 
Fig. [TJd identify the locus where S changes sign. The 
jh = curves have no direct interpretation in terms of 
standard transport coefficients. Finally, if one connects a 
DNLS chain with any two points in the (li, T) plane, its 
profile would correspond to the only path that is charac- 
terized by a constant ratio of j a /jh- 

It is instructive to compare these results with the sce- 
nario expected in the "harmonic" limit, where the nonlin- 
ear terms in the DNLS are negligible. Here, the dynam- 
ics is characterized by an ensemble of freely propagating 
waves and transport is thus ballistic. A direct recon- 
struction of the zero-flux lines by direct simulations is 
not very useful, as, in analogy with the known behaviour 
for the harmonic chain [2!|, the profiles of T and li are 
flat (except for a few sites close to the boundaries). Thus, 
the curves degenerate to single points and no compari- 
son is possible. We thus resort to a different method 
of computing transport coefficients for ballistic systems, 
which is completely analogous to the well-know Landauer 
theory of electronic transport [3(| ■ Consider an iV-site 
chain in between two leads at different temperatures and 
chemical potentials (Tl, lll), {Tr, lir). Since transport is 
ballistic, energy and norm are carried by N independent 
phonon modes, whose dispersion law is U)(q) = 2cosg, q 
being the wavenumber Qq\ < tt). Accordingly, the fluxes 
are iV-independent and the ensuing transport coefficients 
are proportional to N. In this context, the norm and en- 
ergy fluxes are given (up to some numerical constant) by 
the formulae 

r+2 

Ja = dL>t(u>)[f L (u) - f R (U>)\ 



J h = I dwujt(ui)[f L (uj) - f R (uj)] , 

-2 



V. ZERO FLUX CURVES 

A compact pictorial representation of transport prop- 
erties is obtained by drawing the lines corresponding to 



where t(ui) denotes the transmission coefficient, while 
fL,R ar e the equilibrium distribution functions of the 
reservoirs. If we assume that they are composed of 
two infinite linear chains (both with the same disper- 
sion), the equipartition principle implies that the distri- 
butions are of the Rayleigh- Jeans form [l2j], Jl,r(oj) = 
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T 




FIG. 7: (Color online) Zero-flux curves in the (a, h) and (fj,, T) 
planes (panels (a) and (b), respectively). Black dashed and 
blue solid lines correspond to j a = and jh = respectively. 
Simulation are for a chain of length N = 500. The thin dashed 
lines in the upper panel are the T = and T — oo isothermals. 
The thick dot-dashed line identify the locus where S changes 
sign (see text). 



f{TL,R,V>L,R>u) where 



f(T,fJ,,u) 



2(i 



(the factor 2 stems from the definition of z n and Eq. (|3J)) 
The physical meaning of the formulae is pretty intuitive: 
they can be derived from suitable generalized Langevin 
equations 28] following similar steps as for coupled os- 
cillators, see e.g. Ref. [31| . The relevant information is 
contained in the transmission coefficient, that depends on 
how the chain is coupled to the external leads. For the 
Monte-Carlo bath we have used throughout this paper, 
the precise form of t is not known. We thus postulate the 
simplest possible form, namely that, for large TV, t(uS) = t 
for | a; | < 2 and zero otherwise. For our purposes, we set 
t = 1 in the following, otherwise all the coefficients must 
be multiplied by t. If we introduce the function 



*(r,At) 



+ 2 



duj f(T, fj,,oj) = Tin 



which for /! < 
write, 



-1 and T > is always positive, we can 



J a = <f(T L ,^ L )~^(T R ,fi R ) 

Jh = 4(T L -T R )+2(i L $(T L , f i L )-2fi R $(T R , f i R ). 
By expanding to first order in AT = (Tj, — T R ) and 



Afl = (flL 



Mil A// 



where 



J a = M o AT 
J h = Aho AT 

( 4>(7» 
i 4 + 2 M $(T,/i) J^ + 2T$(7» 



(14) 
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With the help of the explicit formulae (fl4j) we can recon- 
struct the zero-flux curves as follows. Starting from an 
initial point (Tj n , /ij„) we compute AT and Afj, inverting 
Eqs. ([M]) setting J a — 0,J h = 1 and J a = l,J h = 0, 
respectively (the value 1 is arbitrary). We then let 
(Ti n , Mm) (Tin + AT, /ij n + A/i) and iterate the proce- 
dure until the whole lines are reconstructed. 

The results are depicted in Fig. HI The curves for the 
linear case are defined only in the region /x < — 1. The 
results of the simulations of the DNLS (solid lines) nicely 
approach the curves of the linear case (dashed lines) upon 
decreasing fi. The agreement is satisfactory, expecially in 
view of the many ad hoc assumptions made in deriving 
Eqs. (O. 



1 ~" — i — 1 — i — ' — i — ' — i — 1 — r 

\ 
\\ 
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FIG. 8: (Color online) Comparison of the zero-flux lines ob- 
tained from simulation of the DNLS equations (black solid 
and blue dot-dashed lines correspond to j a ~ and jh = 
respectively). The dashed (red) lines are the zero- flux lines 
computed by the Landauer formulae as described in the text. 



VI. CONCLUSIONS 

We have presented the first study of stationary trans- 
port properties of the DNLS equation. Due to the non- 
standard form of its Hamiltonian, several new issues have 
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been brought to the fore that had not been addressed 
before in the literature dealing with energy transport in 
oscillator chains. In particular, we have extended the mi- 
croscopic definition of the temperature to the chemical 
potential and defined suitable Monte Carlo thermostat- 
ing schemes to characterize the nonequilibrium steady 
states of the DNLS in various regimes. The simulations 
confirm the expectations that transport is normal (i.e. 
the Onsager coefficients are finite in the thermodynamic 
limit), although some almost ballistic transport is found 
at very low temperature, where the DNLS approaches an 
almost intcgrablc limit. 

Due to the very existence of two naturally coupled 
transport processes, the nonequilibrium steady state can 
display nonmonotonous energy and density profiles. To 
our knowledge, this unusual feature has not been ob- 
served so far in any other oscillator or particle model. 
As seen from Eq. ([TBI) , it is clear that the temperature 
profile cannot in general be linear in y, since the elements 
of A -1 depend on [i and T. In principle, the profiles may 
have have nontrivial shapes depending on the qualitative 
behaviour of the solutions of Eq. (IT51) . In the DNLS, 
the phenomenon is particularly pronounced (the temper- 
ature inside the chain reaches values that are almost three 
times larger than those imposed by the thermal baths) 
because of the strong variability of the Onsager coeffi- 
cients. It would be interesting to find the physical moti- 
vation for this effect to predict and possibly control the 
conditions for its appearance. 

Another novel feature is the fact that the Seebeck co- 
efficient changes sign upon changing the state parame- 



ters e.g. by increasing the interaction strength. The ob- 
servable consequence of this is that the temperature and 
chemical potential gradients change their relative signs. 
As the particle density a increases with /x this also implies 
that a may be larger in the colder regions. 

Furthermore, a remakable feature of the DNLS ther- 
modynamics is the possibility of ne gati ve temperatures 
states in suitable parameter regions [id ]. These regions, 
that are characterized by the presence of long-lived local- 
ized excitations (discrete breathers), have not be consid- 
ered in the present paper, but are definitely worth being 
explored. It may be indeed speculated that they would 
lead to genuine nonlinear transport features and even to 
the birth of new dynamical regimes possibly displaying 
transitions between conducting and insulating states. 

Besides its intrinsic theoretical interest as a testbed for 
the characterization of coupled irreversible processes, the 
DNLS equation opens the way also to experimental in- 
vestigations. In fact, despite its mathematical simplicity, 
the DNLS model can be of guidance in the design and in- 
terpretation of experiments on coupled-transport in cold 
atomic gases in deep optical lattices as well as in optical 
multilayered and nonlinear structures. 
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